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Abstract : 

The functions of the gluon fragmentation into 3 Pj quarkonium are calculated to order 
a 2 s . With the recent progress in analysing quarkonium systems in QCD we show explicitly 
how the socalled divergence in the limit of the zero-binding energy, which is related to P- 
wave quarkonia, is treated correctly in the case of fragmentation functions. The obtained 
fragmentation functions satisfy explicitly at the order of a 2 s the Altarelli-Parisi equation and 
when z — > they behave expected. Some comments on the previous results are 

made. 



1. Introduction 



Quantum Chromodynamics(QCD) is believed to be the correct theory for the strong 
interaction, but its predictive power is limited by the fact that the fundamental quanta, 
the parton, i.e. the quarks and gluons, are not observed freely in Nature. We observe 
in the real world only hadrons. This fact makes it impossible to calculate directly from 
QCD for example an inclusive cross section of single hadron production in a reaction. 
However, in some kinematical regions, the dominant contribution to the cross section can 
be decomposed within the framework of QCD into three parts: the partonic part, the part 
of parton fragmentation into the produced hadron and the part of the parton distributions 
in the initial hadrons. The partonic part can be calculated perturbatively thanks to the 
asymptotic freedom of QCD, while all the long distance physics due to the bound state 
nature of the hadrons is lumped into the parton fragmentation functions and the parton 
distributions. This is the content of the factorization theorem[see Ref. [1] and references 
cited therein]. 

It is hard to study fragmentation functions directly from QCD because of their non- 
perturbative nature. However, the situation changes if the hadron into which partons 
fragment is a quarkonium. It is well known that the bound state effect of a quarkonium 
system can be well described by a nonrelativistic wavefunction if one takes the quarkonium 
system as a color-singlet bound state of a heavy quark Q and its antiquark Q } since in 
the quarkonium rest frame the quark and the antiquark move with a small velocity v. 
Assuming the long distance effect is only contained in the wavefunction, one can calculate 
fragmentation functions perturbatively in QCD[2-8] except a nonperturbative constant: 
the wavefunction at the origin or its derivative at the origin, where one makes an expan- 
sion in the small parameter v. This gives us an opportunity to get an insight into the 
fragmentation functions by starting directly from QCD, hence such study is important. 
The importance can also be shown by the fact that the heavy quark fragmentation func- 
tion for B c meson obtained in this way has similar behaviour as that from the successful 
model for the heavy quark fragmentation[9]. But if one calculates the functions of the 
gluon fragmentation into P-wave quarkonia, the result is infrared (I.R.) divergent at the 
leading order of a s [6,7]. Such I.R. divergence is actually well known in the decay of the 
P-wave quarkonium as the divergence in the limit of the zero-binding energy [10]. The 
appearance of the I.R. divergence indicates that the wavefunction can not contain all the 
long distance effect. It should be kept in mind that the calculations with the mentioned 
assumption are merely calculations within a model, namely the color-singlet model. They 
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are not calculations within the framework of QCD. 

In the color-singlet model one only takes for P-wave quarkonium production the P- 
wave state of the color singlet QQ into account. The probability for production through 
this way is proportional to v 2 . However in [11,12] one realized that a S"-wave state of 
a color octet QQ will contribute to the production rates at the same order of v as the 
P-wave state of the color singlet QQ } because the probability for generating a S"-wave 
state of a color octet QQ is proportional to v° and this state has a probability at order of 
v 2 to be transmitted into a P-wave state of a color singlet QQ through emission of soft 
gluons. This state can also be combined with a soft gluon to form with a probability at 
order of v 2 one of the components of a quarkonium state. From the point of view of a 
relativistic quantum field theory a quarkonium state consists not only of a color singlet 
QQ pair but also of many other components. The S"-wave state of a color octet QQ 
combined with a gluon is a possible component. Therefore one should take the color octet 
QQ into account. Recently, quarkonium systems have been treated rigorously in QCD[12] 
and new factorization formulas for production and decay rates are obtained. Further, 
the nonperturbative effect is represented by matrix elements defined in nonrelativistic 
QCD(NRQCD), not in terms of the wavefunction as in the color-singlet model. For P- 
wave quarkonia there are two nonperturbative constants at the leading order of v instead 
one in the color-singlet model. The second constant is for the contribution from a color 
octet QQ. 

In this work we will calculate with the new factorization formula the functions of the 
gluon fragmentation into 3 Pj quarkonia to order of a 2 , where J = 0, 1,2 is the quantum 
number of the total angular momentum. We assume that the polarization of quarkonia is 
not observed and we will work only at the leading order in v. We will start from the defini- 
tion[13] of the fragmentation functions as in [6]. In Sect. 2 we introduce the definition and 
the new factorization formula for 3 Pj quarkonium. The renormalization of the fragmen- 
tation functions is briefly discussed. The result for the gluon fragmentation at the leading 
order of a s is obtained. For details of the main content in Sect. 2 we refer to Ref[13,12]. 
In Sect. 3 a perturbative calculation in NRQCD for matrix elements is presented. The 
results are needed to extract the perturbative parts in the gluon fragmentation and they 
are universal, i.e. can be used for any process where a P-wave quarkonium is produced. 
In Sect. 4 and Sect. 5 we extend our calculations to order a 2 . Sect. 6 presents our discussion 
and the summary. 

Throughout our work we will use dimensional regularization to regularize ultravi- 
olet(U.V.) divergences and also infrared(I.R.) divergences. With this regularization we 
maintain the gauge-invariance. 
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2. The Definition of The Gluon Fragmentation Function 
and The Leading Order Result 

We give the definition of the gluon fragmentation function in d dimensions. To give the 
definition it is convenient to work in the light-cone coordinate system. In this coordinate 
system a (/-vector p is expressed as p^ = Pt), with p + = (p° + p d_1 )/\/2, p~ = 

(p° — p d ~ 1 ) I \/2. Introducing a vector n with n' 4 = (0, 1, 0, • • • , 0) = (0, 1, Ot), the gluon 
fragmentation function for a spinless hadron H or for a hadron without observing its 
polarization is defined as[13]: 

n(°) ( z ) - - / dr~e~ iP x lz < Q\G b ' +v (Q} 

r-OO 

{Pexp{-ig s / dXn ■ G(\n»)}} hc cJ H (P + ,0 T )a H (P + ,0t) (2-1) 
Jo 

r-OO 

{Pexp{ig s / d\n ■ G{\n»)}} cd G d >t{n,x- ^ T )\<d >, 

J x~ 

where G^(x) = G^(x)T a , G a (x) is the gluon field and T a (a = 1, ... ,8) are the color ma- 
trices. G a,ilv is the gluon field strength tensor and aJj(P) is the creation operator for the 
hadron H . For hadrons with nonzero spin the summation over the spin is understood. The 
definition is a unrenormalized version. In calculating D^ G (z) one will encounter U.V. di- 
vergences, requiring renormalization. Following [13] the renormalized gluon fragmentation 
function can be defined as: 




Here the function L a (z) is chosen so as to cancel the U.V. divergences. In MS scheme 
L a (z) takes the form: 

L a (z) = Y,-^Li N \z), (2.3) 

N 

where e = 4 — d. From Eq.(2.2) one can derive the Altarelli-Parisi type of the evolution 
equation for the fragmentation function. We will use the modified MS scheme. That means 
in practical calculations that the terms with N e = -| — 7 + ln(47r) will be cancelled by the 
function L a (z). The function Dh/q(z) is interpreted as the probability of a gluon G with 
momentum k to decay into the hadron H with momentum component P + = zk + , it is 
gauge invariant from the definition. Further the function is also invariant under a Lorentz 
boost along the moving direction of the hadron and under a rotation with the direction as 
the rotate axis. 
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Now we give the new factorization formula for 3 Pj quarkonium production. We use 
the notation \J f° r the 3 Pj quarkonium. A production rate for \J takes the following 
factorized form [12]: 

<XJ) = < 0\OrCPj)\0 > +^0^ < OIOFCSJIO > . (2.4) 

The operators 0* J ( 3 Pj) and Og J ( 3 Pi) are given by: 
Of ( 3 Si) = x t ^T>(a| f a Jf )V' t ^T a X 

o?( 3 Po) = \xH-^v ■^(«W^ t (- ^ 6 .*) x 

o? ( 3 po = - ^ 5 x,).^(aW^(- i 5 x,) iX (2 ' 5) 

Of ( 3 P 2 ) = xH - \ D {1 a j} )^(a^a H )^( - % - D {l a j} ) x , 

where D is the space part of the covariant derivative D^ 1 and a^i = 1,2,3) is the Pauli 
matrix. T a (a = 1,...,8) is the SU(3) color matrix. The notation {ij} means only to 
take the symmetric and traceless part of a tensor. In Eq. (2.5) tp and x^ are fields with 
two components for the quark Q and the antiquark Q in NRQCD. M is the mass of the 
quark. a~j^ is the creation operator for the hadron in its rest frame. The matrix elements in 
Eq.(2.4) are defined in NRQCD. In Eq.(2.4) the part with Pg is the contribution from the 
color octet QQ in a 3 Si state. The matrix elements represent the nonperturbative effect, 
while the coefficients Pi and Pg can be calculated perturbatively and they should be free 
from I.R. singularities. 

For calculating the coefficients Pi and Pg one can takes an on-shell QQ pair with 
small relative momentum instead of the hadron in Eq. (2.4), and calculates a(QQ) with 
perturbation theory. Then one also calculates the matrix elements in NRQCD with the 
same QQ pair. Substituting the results into Eq.(2.4) one can identify the coefficients. This 
is called matching. There will be I.R. divergences in a(QQ) and also these I.R. divergences 
will appear in the matrix elements so that the coefficients are free from them. 

According to Eq.(2.4) we can write the form for the gluon fragmentation function as: 

D Xj /g(z) = D^z, J)^ < 0\OrCPj)\0 > +Ds(z)^ < 0|O 8 ^( 3 Si)|0 > 
= D Xj/ga (z) + D Xj/Gj8 (z). 

The function Di(z,J) and D%{z) are dimensionless. In Eq.(2.6) the contribution from a 
color octet QQ pair and from a color singlet QQ pair is given by D XJ /g,s( z ) an d D XJ /g,i( z ) 
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respectively, with corresponding perturbative parts D$(z) and Di(z } J). We may call them 
the color-octet and color-singlet components. The function D$(z) is same for J = 0,1,2. 
The task of the present work is to determine D%{z) and £>i(z, J) to order of a 2 s . 

At order of a S} only D$(z) is nonzero, while Di(z) becomes nonzero at order of 
a 2 s . At leading order of a s there is one Feynman diagram in Fig.l which contributes to 
D XJ /g,9,{ z )- I n Fig. 1 the double line presents the line operator in Eq.(2.1), the broken 
line is the Cutkosky cut. The calculation can be easily done with the matching procedure 
mentioned before. We obtain 

Ds(z) = ^^m-^) (2-7) 

This result is independent of d. To calculate contributions in the next order of a S} we 
need to calculate the matrix element < 0|O^( 3 Si)|0 > to order of a s , where the hadron 
is replaced by a QQ pair as required by the matching. We will do the calculation in the 
next section. 



3. The Matrix Element in NRQCD 



At the tree-level, the contribution to the matrix element < 0|Og ( 3 Si)|0 > is pre- 
sented by the diagram in Fiq.2a, where the black circle means the vertex sandwiched 
between x an d m the operator. For the matrix element here it is <JiT a . The line means 
the quark line. The broken line is the Cutkosky cut. We let the QQ pair on the left side 
of the diagram move with a relative velocity 2v and the QQ pair in the right side of the 
diagram move with a relative velocity 2v'. The QQ pair is in CMS. At the tree-level we 
obtain 

< OlOf^SOlO >= rfaiT'&aiT'ri, (3.1) 

where rj^ and £ are the Pauli spinors with two components in NRQCD for the quark and 
antiquark respectively. 

In the next order there are two types of the diagrams contributing to the matrix 
elements. One type is given in Fig. 2b. This type of the diagram includes the correction to 
the vertex and the correction to the external quark line. The other type is given in Fig. 2c, 
where a real gluon is emitted by the quarks. This causes the transition from a color octet 
QQ into a color singlet QQ. 

Using NRQCD we calculate the contribution from these two types of the diagrams. In 
the calculation we make an expansion in v = |v| and v' = |v'| and neglect in the expansion 
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irrelevant higher order terms. The contribution from the type of the diagrams in Fig. 2b 
before the renormalization is: 



M Fig .2 b = j^JaiT'tfviT'ri ■ {(3 - ^(v 2 + v' 2 ){±- + ± + ln(-^)) 



7T 2 , , IT 2 



(3.2) 



( 1 +- 2 )-^7( 1 + -' 2 )} 



2Av v ' 2Av 



where ^ is the U.V. divergence and \i is the corresponding scale, jj means the I.R. diver- 
gence with ej = d — 4 and /ij is the corresponding scale. We will keep this notation in our 
work. The terms with - and -\ are the Coulomb singularities. 

V V ° 

The contribution from the type of the diagrams in Fig. 2c is: 



MFi g .2c = 7 ^(- + -+H JL )) 

(2tt 2 ) v ej e /ij 



• { - (3 + (v 2 + /j^a.r^V.r,] + • v '( Va^ 1 ^ ( 3 - 3 ) 
+ ^v i T a tfv i T a ri)}. 

In Eq.(3.3) the terms with v • v' in the second line can be identified as the matrix elements 
of the operators Of®( 3 Pj). These terms present the transition mentioned before. For 
renormalization one introduces for the operators a matrix of the renormalization constants 
because of the operator mixing. After renormalization the result is free from U.V. diver- 
gence. The I.R. divergence still remains and it represents the nonperturbative nature of 
the matrix element. We obtain: 

< o|of^)|o > = {i - ^(^ + ^)} < 0103^(^010 >o 

+ j£y ■ ^ • ^ < o|o?*( 3 p ) + o?*( 3 po + o?*Cp 2 )\o >o 

— + -(7 - ln(4vr) - ln(-^-) + ln(-^)l + • • • , 

(3.4) 

where the • • • is for terms in higher orders in g S} in the velocity v and v' and also for 
a term related to a operator which is irrelevant in this work. The subcript o with the 
matrix elements means only to take their tree-level results. A by-product generated in the 
calculation is the renormalization group equation for the matrix element < 0|Og J ( 3 Si)|0 >• 
It reads: 

/<°' ^"°> ^.(„^<0| »(^)|0> ( 3,) 
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This equation was first derived in [12,11]. 

The result in Eq.(3.4) is useful not only for our calculation in the following sections 
but also for 3 Pj production which may be not due to fragmentation. Using this result we 
can now calculate higher order contributions to the fragmentation functions in Eq.(2.6). 
At higher order of a s it may be difficult to perform the matching procedure mentioned 
in Sect. 2. It is still convenient to use wavefunctions to project out different states from 
a QQ pair. However, the final results for L>i(z,J) and D%{z) do not depend on these 
wavefunction. Details about such projections may be found in [14]. We introduce a 
radial wavefunction Ri(r) to project the 3 Pj color singlet QQ state and a octet radial 
wavefunction R^ a \r)(a = 1,...,8) to project the 3 S\ color octet QQ state. With these 
wavefunctions we obtain for the matrix elements at the tree-level in Eq.(2.5): 

< o\orcpj)\o > = 9(2 ^ +1 V i(o)i 2 , 

<o\orCs 1 )\o> = ^- £l4 c) (o)l 2 , 

c 

where R'i(0) is the first derivative of i?i(r) at the origin. The result in Eq.(3.4) can be 
expressed with these wavefunctions as: 

< oio^soio > = ^J2 l4 c) (o)| 2 (i - < >) 

c ^ ' 

I 9 2 s 16 1 9( 2 ^ + V (0) | 2 (3.7) 
+ (2vr)2 27 M? 2vr |jKl(Ujl 

I — + -(7 - ln(4vr) - ln(-^-) + ln(-^)|. 
The quantity < u _1 > is the average of u _1 with the color octet wavefunction. 



4. The Color-Singlet Component 



In this section we explain our calculation in detail for J = 1 and simply give the 
results for J = 0,2. As mentioned in Sect. 3 we will use the color singlet wavefunction to 
project out the 3 Pi color singlet QQ state. In order to present the Feynman diagrams in 
the calculation efficiently, we decompose the function Dh/q(z) as: 

D hig {z)~T h ®T]j. (4.1) 

Here the symbol ® denotes contractions of color- and Lorentz-indices and summations over 
all possible intermediate states. Such a decomposition is always possible by sandwiching 
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the operator 1-^" X -^1 between ajj and ajj in the definition in Eq.(2.1). Tjj may be 
called fragmentation amplitude. 

The color-singlet component D Xj /qi (z) becomes nonzero at the order of a 2 . The 
Feynman diagrams for Tjj are given in Fig. 3. Since the QQ pair is in a color singlet state, 
there must at least be two gluon lines attaching the quark line. There is a I.R. divergence 
if the energy of the gluon exchanged between two quarks becomes zero. This results a 
singularity in D XJ /g,i( z ) at z = 1. We use dimensional regularization to regularize this 
I.R. divergence. The calculation is tedious but straightforward. We obtain: 



7T 2 



• I (- + -( 7 - ln(4vr)) - ln(-^) + - Wl - z) + — L- - 1 - -z - z 2 } . 

(4.2) 

Here the I.R. singularity comes from the term (1 — z)~ 1+ei . We used the expansion in ej 
for the distribution (1 — z)~ 1+ei 

,1 - = > - + ^ + + • • • , («) 

and neglected termes in higher order in ej. The + prescription is the same as in [15]. 
To extract the perturbative part L>i(z,J = 1) in Eq. (2.6) one notes that there is a 
contribution from the color singlet wavefunction in < 0|Og 1 ( 3 Si)|0 > at the one-loop level 
in Eq.(3.7) and we obtain at the order of a 2 s or g^: 

27., 2 1 



D Xl/GA (z) = D 1 (z,J = l) — \R' 1 (0)\ 



2tt' iV n M 5 

+ ^("'w-i-^ |fl '' (0)|2w -^ < 4 ' 41 

{ ± + » (7 _ M 4,)- M |L ) + M ^ )} . 

Comparing Eq.(4.4) and Eq.(4.2) we note the I.R. divergence in the r.h.s. of Eq.(4.4) 
appears with the same coefficient in the l.h.s., i.e. in Eq.(4.2). The physical meaning here 
is clear. The QQ pair before emission of the real gluon is in a S"-wave color-octet state. 
When the gluon becomes soft, the color octet QQ pair becomes real and the emission of 
the soft gluon is nonperturbative. This nonperturbative effect is correctly represented by 
the matrix element < 0|Og 1 ( 3 Si)|0 >, and hence D\(z, 1) is I.R. finite and it is: 

Di(z, 1) = —OL 2 Aa) i [- - ln(-^)U(l - z) + 1 —— - 1 - -z - z 2 \ . (4.5) 
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Performing the cacluation for J = 0, 2 we also obtain D\(z, J = 0, 2), which are finite, 

d^,j = o) = [\ - M^)W - + 

+ -(5 - 32)ln(l -z) + -(-26z 2 + 85z - 8)} 

4 8 (4.6) 
9 19 // 1 v ; 

V ' ; 81 sWU 40 K 2M n K J (l-z) + 

+ ^(2 - z)ln(l -z) + \(-±z 2 + \\z - 4)}. 

5 4 

The results in Eq.(4.5) and Eq.(4.7) are actually not new. They were first obtained in 
[6] for J = 1 and in [7] for all J. However, the results here do not agree with those in 
[6,7] at z = 1. For z ^ 1 they agree. The disagreement clearly depends on how the I.R. 
divergence is regularized and some comments are in order. In [6] a cutoff for the transversal 
momentum of the exchanged gluon is introduced for the regularization and the singularity 
is absorbed by the matrix element through a renormalization group equation. As we can 
see from the calculation in Sect. 3, such cutoff is meanless since the transversal direction can 
not be defined and the introduction of the cutoff violates the rotation invariance. Therefore 
the result in [6] is not correct for z = 1. In [7] a cutoff A of the gluon energy is introduced, 
and this cutoff is regarded as a scale for separating short distance effects at the energy scale 
M from long distance effects at the scale Mv or smaller in quarkonium systems. Hence in 
the result from [7] there is a term for the I.R. singularity like ~ 6(1 — z)lnA. Such term 
is actually cancelled if one takes for the matrix element the one-loop result, where the 
same cutoff should be used. The effect may be only to replace A with the renormalization 
scale /i. However such a cutoff can violate the gauge invariance and the invariance of a 
Lorentz boost in the hadron moving direction. The disagreement may be due to these 
violations. We used the dimensional regularization and the gauge invariance is respected 
in each step. Generally, extracted perturbative results are different by separating short 
distance effects from long distance effects through introducing explicitly a separating scale, 
which usually is some cutoff, or through identifying I.R. divergences as long distance effects 
in dimensional regularization, since the definition of nonperturbative parts, which although 
is hard to give, is different. An example may be found in theoretical attempts to explain 
EMC effect [16]. In [7] the coefficient for the term with 6(1 - z) is ±§ - In ^, |f - In A 
and — In for J = 0, 1 and 2 respectively. Comparing these with our corresponding 
results | — In f — In and ^| — In one can see that the disagreement is different 
for different J. This means that one can not relate the cutoff A to /i with a simple relation 
to solve the discrepancy and the discrepancy can not be absorbed into the definition of the 
color-octet matrix element. 
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5. The Color-Octet Component 



In this section we proceed to calculate the color-component D XJ /g,8( z ) } or D$(z) 
upto order of a 2 . This component is independent of the total angular momentum of the 
quarkonium. We use the color octet wavefunction to project out the color octet QQ state. 
According to the structure of the Feynman diagrams we write D XJ /g,s( z ) with the result 
in Sect. 2 as the following form: 



;w;) = i*(i-^D^(°)i 2 ^ (M) 

+ D A (z) + D B (z) + Dc(z). 

Here we neglected terms of order higher than a 2 . The term in the first line is the result 
from Sect. 2. The function Da(z) } Dp(z) and Dc{z) will now be calculated 

Da{z) is the contribution by adding the conventional one loop correction for the 
quark-gluon vertex, the gluon propagator and the external quark line in Fig.l. The U.V. 
divergences are simply cancelled by the renormalization constants of the coupling and the 
wavefunctions. We obtain 



B A {z) = -l±-^y |i?8(0)| 2 — rr— - • { (- - -N f ) In ( 



(4^)28^^' 8V yi M 3 16 lv 3 9 " HM 2 
-^<„-.> + (l + 7 - N 4.,-h 1 (J r )) 

Nf 9 

A J Tfl 

+ 1.822+ -^ReM^)}^!-,), 



(5.2) 



Fp(x) = / dot(ot — a 2 ) In (a 2 — a + x — i0 + ^ 
Jo 



where Nf is the total number of the quark flavours in QCD. The Coulomb singularity 
< u _1 > comes only from the vertex correction. The term with the function Fp comes 
from the correction of the gluon propagator through quark loops. This function will have 
an imaginary part if x < 1 and only the real part contributes. The number 1.822.. comes 
from some factors like In 2 and In 3. In Da(z) there is a I.R. singularity. 

Db{z) is the contribution from the Feynman diagrams given in Fig. 4a for the frag- 
mentation amplitude Th ■ This contribution is proportional to 6(1 — z). In this part there 
is a U.V. divergence and also I.R. divergences as the poles e~ 2 and e _1 . Here some details 
in the calculation may be needed to give. In the calculation we have such type of the 
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integral over a loop momentum q 

( d( l ,d 1 1 1 [ r , n . fl( U ■ g) ) /,o^ 

l 2vr J (P + g) 2 + z0+ ' q 2 + z0+ ' -n ■ q + z0+ ' 9J P • q + g 2 + z0+ j 1 ' J 

where P has only + and — components and P 2 = 4M 2 . The term with fo is U.V. divergent, 
whereas both terms are I.R. divergent. To calculate the integral one can perform first the 
integration over q~ via Cauchy's theorem and then integrate over qx- The U.V. divergence 
will appear after the qx integration as r(|). In order to identify various divergences we will 
keep the factor (d — 2) and factors like T-function combined with 7r, which are generated by 
the qx integration. The remained factors are expanded in e. The result for the contribution 
before renormalization is: 

n (o) r ^_ 9$ 3 V^ipM/n^ 1 1 2 (a-\^( V \ A ~ d 

+ (-7 + i + 4 + y)- r (|)W 1 -^- 

Here we neglected irrelevant terms in higher order of e. The singularity in T(|) is the U.V. 
one. Other singularities are I.R. singularities. 

Now we turn to the last part Dc{z). This part receives contributions from the dia- 
grams in Fig. 3 and also from those in Fig. 4b. In this part there is a real gluon exchange 
and we have to integrate over its transversal momentum. The integration over the + com- 
ponent and — component of the gluon momentum is easily done by the (^-functions, one of 
which is the on-shell condition and the other is due to the total momentum conservation 
in the + direction. Again there are U.V. divergence and I.R. divergences. We obtain: 

° c {Z) - (I^2J^ M^16^2 (47r) W 

• { (4*(i - *) - *(i - *) + ^r^ 1 )^ + - - 2 - 2 + - 3 )) • r(i + f ) 

6 L ~ z 2 (5.5) 

+ ( - - + 26(1 -z)+ _ z'^l -2z + 3z 2 - 2z 3 + z 4 ) 
6 z )+ 

ijY^ + (l^)J^- 1 (i - 2z + 3, 2 - 2, 3 + z-)) • r(|)}. 



The U.V. divergence is r(|). Other poles in e are I.R. singularities. Adding D^\z) 
and D^j\z) together the I.R. singularity with e -2 is compensated. Further, the I.R. 
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singularity in the coefficients in the front of is cancelled too. Now we can safely 

perform renormalization according to Eq. (2.2) with the result we obtain 



+ Dc(z) = ^ylY. l4 C) (0)| 2 ^ • { [S(l - z) 

+ 2 z-\\ ~2z + 3z 2 - 2z 3 + z 4 )] In (- ^ 



(1-*)+ J v 4M 2y 

2 u 2 ( 5 - 6 ) 

+ [-- + ln^- 1 + ln(J^)]6(l-z) 



2vr 2 

+ (41n2 + — )8(1 -z) 



- z~ 1 {2 -Az + z 2 - z 3 )ln(l - z)} 

where we still have a I.R. singularity ej 1 . However this singularity will be cancelled if we 
add Da{z). The final result for D XJ /g,s(z) is: 

J\, /0 ..W = j|«(i-*)£E i^W^ 

c 

+ ^.-'(1-2, + 3.-' -2.-' + .-«)] In (^L) <") 



- z _1 (2 - 4z + - z 3 )ln(l - *) 

+ [ - - < > +4.4946 + — + 3 E R e ^(^)] *(1 " *)} • 



To extract D%{z) we use the result of the quark matrix element in Sect. 2 and obtain the 
equation for D%(z) to order a 2 

(5.8) 



+ (M')-&ta-*))±'Z\rt\o)\*-± 



96 v /; 87r ^ 1 8 Wl M 3 ' 

c 

We note that the Coulomb singularity in Eq.(5.8) appears with the same coefficient as in 
Eq.(5.7). This fact results that D$(z) is free from this singularity. From Eq.(5.8) and (5.7) 
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one finds that 



Ds(z) = ^a s (n)S(l - z) 

16 uv 3 9 fJ K ' 



+ jY^7j; z ~ 1{1 ~ 2z + 3z2 ~ 2z3 + z4)] ln (iM^) ( 5 - 9 ) 

- z~ 1 {2 -4z + z 2 - z 3 )ln(l - z) 

2 N f 2 

+ [4-5946 + \ + 1 E Re M^)]^(l " *)}• 

Now we complete our calculation for gluon fragmentation into 3 Pj-quarkonia. The frag- 
mentation function is 

D xj/G (z, = Mz, J) < 0\OrCPj)\0 > ±r + D s (z) < OWCSJIO > ± (5.10) 

with D\(z, J) and D%{z) given in Eq.(4.5), Eq.(4.6) and Eq.(5.9). Each component is now 
known to order a 2 . The nonperturbative effect is included in the two matrix elements. 
It is interesting to check whether the fragmentation functions calculated here satisfy their 
evolution equation. For gluon fragmentation the evolution equation reads 



^ h/g{ = f < ^p G ^ G ( z / y ^^ DH/G (y^^ + ^2 I —P G ^ q (z/y } fx)D H/q (y } fi) 

*y z y q *j z y 

(5.11) 

In one loop approximation the splitting functions Pq^g(x } //) and Po^ q (x } /i) are the same 
as these in the evolution of parton distributions: 

P G ^g(x, Ll) = < : 1 1 + X - X 2 + ( J -)S(1 - z) } , 

V ' HJ vr \(l-x) + x V 12 18 ; V ; J' 



P G ^ q (x^) = ^{x 2 + (l-x) 2 ). 
Taking the derivative of our gluon fragmentation function we get: 



(5.12) 



d/i 



-(z,») =< 0|O-( 3 ^)|0 > • { [ y - V8 N f ]S(l - z) 

+ , _ 2 . z-\l -2z + 3z 2 - 2z 3 + z 4 )}. 
I-*- z )+ 



(5.13) 



The //-dependence in the color-singlet component L>i(z,J) is compensated by that in 
the matrix element of Og J ( 3 S\) through the renormalization group equation in Eq.(3.5). 
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At order of a 2 the nonzero contribution to the r.h.s. of Eq.(5.11) is only from gluon 
fragmentation function because quarks can not fragment into the quarkonium at order 
of a s . Substituting the lowest order result of the gluon fragmentation function into the 
r.h.s. in Eq.(5.11) we find after a little algebra, that our fragmentation function satisfies 
Eq.(5.11) at order of a 2 . Here we make some comments on the results for the fragmentation 
of c-quark into \cJ given in [8]. Since a gluon can decay into 3 Pj quarkonia at order of 
a s , there must be a //-dependence in the function of the quark fragmentation into 3 Pj 
quarkonia at order of a 2 according to the Altarelli-Parisi equation. However, from the 
results in [8], there is no //-dependence at this order. Hence the results in [8] for c quark 
fragmentation into XcJ are not correct. A new calculation is needed and is already on the 
way [17], where it is realized that light quarks can also decay into 3 Pj quarkonia through 
the color octet QQ state at order of a 2 . 

From our results we can see that the gluon fragmentation functions behave like z~ x 
as z — > 0. This is expected in general or simply from the evolution equation (5.11) and 
(5.12). However, our results tell us that at the scale /i = 2M such behaviour disappears. 

6. Discussion and Summary 

To discuss possible impact of our results for 3 Pj quarkonium production, we calculate 
the first moment of the fragmentation functions. This moment is the probability to 

find a quarkonium without observing its momentum in a jet initialized by a gluon. A rough 
estimate for 3 Pj quarkonium production due to fragmentation can be made by taking the 
product of parton cross sections and the first moments of the partons. We take the c 
quark as an example. We neglect the heavy flavor, i.e. b and t quark and take Nf = 4 and 
m q = for light flavors. To avoid the singularity at z = 0, we calculate the moments at 
fj, = 2M C . We obtain: 

M<°>(G- Xc0 ) 



M<°>(G- Xc i) 
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M 3 



< 0|Of°( 3 5i)|0 > -(0.131a a (2M c ) + 0.537a^(2M c )) 



1 



0.067^(2M C )— < 0|Of °( 3 P )|0 >, 

c 

1 



M 3 



< 0|Of 1 ( 3 5i)|0 > -(0.131a a (2M c ) + 0.537a^(2M c )) 



1 



0.021^(2M C )— < 0|Of 1 ( 3 P 1 )|0 >, 



M (0) (G ^ Xc2 ) = Jp< O|0f 2 ( 3 Si)|0 > -(0.131a a (2M c ) + 0.537^(2M C )) 

(6.1) 

- 0.043a*(2M c ) — < O|0f 2 ( 3 Pi)|0 > . 

It is interesting to note that the color singlet component gives a negative contribution. 
The matrix elements are until now not known precisely. For different J they are related to 
each other through a spin factor if one neglects higher order terms in v in NRQCD action. 
A rough estimation for these matrix elements can be found in [11,12]. Using the results 
there and M c = 1.5GeV, a s (2M c ) = 0.26, the values of the first moments are 0.9 X 10 -4 , 
3.8 X 10 -4 and 5.5 X 10 -4 for \c0, Xa and \c2 respectively. The matrix elements may be 
possibly determined in lattice study. Their corresponding matrix elements in decays have 
been studied on lattice in [18]. However it is problematic to study the matrix elements 
here for production, because in spite of the complicated structure there is a sum of all 
possible intermediate state combined with a quarkonium. It remains to see how such 
matrix elements can be studied on lattice. 

In this work we have shown how the infrared singularity, which appears in the color 
singlet model, is correctly treated in the case of the gluon fragmentation functions in 
the context of recent progress. In the color-singlet model the nonperturbative effect in 
3 Pj quarkonia is not well separated and some important contribution is missing, thus 
the infrared singularity appears. In the rigorous treatment in the framework of QCD, 
such singularity disappears as we have shown. Further, the Coulomb singularity is also 
absorbed by the matrix elements. The final results we obtained are free from infrared 
singularities and Coulomb singularities. A clear separation between long distance effect 
and short distance effect is achieved. In old calculations within the color-singlet model 
for 3 Pj quarkonium production and decay such separation can not be established and 
the contribution from color octet QQ pairs is neglected. In that sence, new calculations 
must be done in the way shown here. The results in Sect. 3 may be needed in these new 
calculations. It is interesting to note a color-octet component can receive contributions 
in a production at lower order of a s than a color-singlet component. This fact may be 
useful for example to explain the discrepancy between theory and experiment for direct xi 
production at hadron colliders[see [19] and references cited therein], where the important 
partonic process, i.e., \i production via gluon fusion, proceeds at order of a 3 in the color- 
singlet model. If one takes the contribution from color octet QQ pairs, the process can 
happen at order of a 2 s . This maybe explains why the production rate predicted within the 
color-singlet model is lower than indicated by experimental data. It should be pointed out 
that the Landau- Yang theorem does not apply here because gluons have color as extra 
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freedom than photons have. 

Although full calculations for production rates of quarkonia can in general be per- 
formed, but in cases where the factorization theorem applies the amount of work can be 
greatly reduced by using the concept of fragmentation. One needs only to calculate parton 
production rates and use fragmentation functions. Nevertheless the importance of our 
work is not limited only to this issue as a good approximation. As already pointed out 
in the introduction, with quarkonia one can start directly from QCD to study properties 
of fragmentation functions and the results can be used to build models for fragmentation 
into hadrons consisting of one heavy quark with light quarks. The gluon fragmentation 
functions are generally expected to be singular as z — > 0. Through our work here, the 
obtained fragmentation functions are regular at the whole range of z at the scale /i = 2M. 

In the literature there are different methods for calculating fragmentation functions. 
In [3,5] a method is proposed by separating a fragmentation amplitude and a parton 
production amplitude in a full amplitude of quarkonium production. The separation is 
achieved at the tree-level by choosing the light cone gauge: n ■ G = 0[5]. This method 
is equivalent to the method here. If we choose the light cone gauge, the path-ordered 
exponential of the gluon field in the definition of a parton fragmentation function is unity, 
and we can obtain the same Feynman diagrams and the same fragmentation amplitudes 
presented in [3,5]. In the light cone gauge the number of Feynman diagrams can be reduced 
considerably. In our case only one diagram in Fig. 4b remains in this gauge while the others 
in Fig. 4a and Fig. 4b do not exist anymore. But it is not convenient to work with this gauge 
because the divergences are very difficult to handle. In our work we have taken Feynman 
gauge. 

To summarize: We have calculated to order a 2 s the functions of gluon fragmentation 
into 3 Pj quarkonia, the results are given in Eq.(4.5), Eq.(4.6), Eq.(5.9) and Eq.(5.10). 
The perturbative parts of the functions are free from infrared and Coulomb singularities. 
These singularities are correctly factorized into a set of matrix elements, which are defined 
in NRQCD and represent nonperturbative effect. The obtained fragmentation functions 
satisfy the Altarelli-Parisi evolution equation explicitly as they must. An useful result for 
a matrix element in NRQCD is also obtained, which may be needed to calculate cross 
sections of single xj production. 
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Figure Caption 



Fig.l: The Feynman diagram for the color octet component at the leading order of a s . 
Fig. 2a: The Feynman diagram for matrix element of the operator 0®® ( 3 Si) at the leading 
order of a s in NRQCD. 

Fig. 2b: Some Feynman diagrams for a part of the contribution to the matrix element at 
the order of a s . In this part only the vertex correction and the external leg correction is 
included. 

Fig. 2c: Some Feynman diagrams for a part of the contribution to the matrix element at 
the order of a s . In this part only the contribution induced by a real gluon exchange is 
included. 

Fig. 3: The Feynman diagrams for the contribution to the color singlet component at the 
order of a 2 s . They also lead to a part of the contribution to Dc(Z) in Eq.(5.1) for the color 
octet component. 

Fig. 4a: The Feynman diagrams for the contribution to Db{z) in Eq.(5.1). 
Fig. 4b: The Feynman diagrams for the contribution to Dc{z) in Eq.(5.1). 
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Fig.2a 



Fig. 2b 



Fig. 2c 




Fig.4a 




Fig.4b 
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